/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

package com.bioinfseq.bioinfmultiseq;

import java.io.BufferedReader;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.logging.Level;
import java.util.logging.Logger;
import org.biojava3.alignment.Alignments;
import org.biojava3.alignment.GuideTree;
import org.biojava3.alignment.SimpleSubstitutionMatrix;
import org.biojava3.alignment.template.PairwiseSequenceAligner;
import org.biojava3.alignment.template.PairwiseSequenceScorer;
import org.biojava3.alignment.template.Profile;
import org.biojava3.core.sequence.DNASequence;
import org.biojava3.core.sequence.compound.AmbiguityDNACompoundSet;
import org.biojava3.core.sequence.compound.DNACompoundSet;
import org.biojava3.core.sequence.compound.NucleotideCompound;

/**
 *
 * @author Bartek
 */
public class App {

    public static void main(String[] args) {
        String inputPath = "seq_input.txt";
        BufferedReader reader = null;
        int length = 0;
        double[] scores;
        boolean withMatrix = false;
        String firstSequence = "", secondSequence = "", option = "", penalties = "";
        short openPenalty = 0, extensionPenalty = 0;
        List<DNASequence> sequences = new ArrayList<>();
        SimpleSubstitutionMatrix<NucleotideCompound> substitutionMatrix = null;
        Profile<DNASequence, NucleotideCompound> profile;
        
        try {
            reader = new BufferedReader(new FileReader(inputPath));
            String sequence;
            while ((sequence = reader.readLine()) != null) {
                if(sequence.equals("Matrix")){
                    withMatrix = true;
                    break;
                }
                if(sequence.length() > length)
                    length = sequence.length();
                sequences.add(new DNASequence(sequence));
            }
            if(withMatrix)
            substitutionMatrix
                        = new SimpleSubstitutionMatrix<>(DNACompoundSet.getDNACompoundSet(), reader, "mm");
        } catch (FileNotFoundException ex) {
            Logger.getLogger(App.class.getName()).log(Level.SEVERE, null, ex);
        } catch (IOException ex) {
            Logger.getLogger(App.class.getName()).log(Level.SEVERE, null, ex);
        } finally {
            try {
                if (reader != null) {
                    reader.close();
                }
            } catch (IOException ex) {
                Logger.getLogger(App.class.getName()).log(Level.SEVERE, null, ex);
            }
        }
        
        CustomPenalty cPenalty = new CustomPenalty(openPenalty, extensionPenalty);
        
        for(int i = 0; i < sequences.size(); i++){
            for(int j = 0; j < sequences.size(); j++){
                PairwiseSequenceAligner<DNASequence, NucleotideCompound> localAligner
                    = Alignments.getPairwiseAligner(sequences.get(i), sequences.get(j),
                            Alignments.PairwiseSequenceAlignerType.LOCAL, cPenalty, substitutionMatrix);

                PairwiseSequenceAligner<DNASequence, NucleotideCompound> globalAligner
                    = Alignments.getPairwiseAligner(sequences.get(i), sequences.get(j),
                            Alignments.PairwiseSequenceAlignerType.GLOBAL, cPenalty, substitutionMatrix);
                
                printAlignment("Local alignment", localAligner);
                printAlignment("Global alignment", globalAligner);
            }
        }
        
        
                
        if(withMatrix)
            profile = Alignments.getMultipleSequenceAlignment(sequences,substitutionMatrix);
        else
            profile = Alignments.getMultipleSequenceAlignment(sequences);                
        
//        scores = Alignments.getAllPairsScores(sequences, Alignments.PairwiseSequenceScorerType.GLOBAL, cPenalty, substitutionMatrix);
//        System.out.println();
//        System.out.println("Pairs Scores:");
//        for(int i = 0; i < scores.length; i++){
//            System.out.print(String.valueOf(scores[i]) + " ");
//        }
        
        List<PairwiseSequenceScorer<DNASequence, NucleotideCompound>> scoresList = 
                Alignments.getAllPairsScorers(sequences, Alignments.PairwiseSequenceScorerType.GLOBAL, cPenalty, substitutionMatrix);
        System.out.println();
        System.out.println("getAllPairsScorers():");
        for(int i = 0; i < scoresList.size(); i++){
            System.out.println(scoresList.get(i).getScore());
        }
        
        Alignments.runPairwiseScorers(scoresList);
        GuideTree<DNASequence, NucleotideCompound> tree = new GuideTree<DNASequence, NucleotideCompound>(sequences, scoresList);
        
        System.out.println();
        System.out.println("GuideTree:");
        System.out.println(tree);
        
        Profile<DNASequence, NucleotideCompound> msa;
        msa = Alignments.getProgressiveAlignment(tree, Alignments.ProfileProfileAlignerType.GLOBAL, cPenalty, substitutionMatrix);
        
        System.out.println();
        System.out.println("getProgressiveAlignment():");
        System.out.println(msa);
        
        System.out.println();
        System.out.println("Multialignment:");
        System.out.println(profile);
        System.out.println();
        System.out.println("Profile:");
        
        DNACompoundSet dnaSet = DNACompoundSet.getDNACompoundSet(); 
        List<NucleotideCompound> compounds = new ArrayList<>();
        compounds.add(dnaSet.getCompoundForString("A"));
        compounds.add(dnaSet.getCompoundForString("G"));
        compounds.add(dnaSet.getCompoundForString("T"));
        compounds.add(dnaSet.getCompoundForString("C"));
        compounds.add(dnaSet.getCompoundForString("-"));
        //List<NucleotideCompound> compounds = dnaSet.getAllCompounds();
        
        System.out.println("   A  G  T  C  -");
        for (int i = 1; i <= length; i++) {
            System.out.println( String.valueOf(i) + ": " + Arrays.toString(profile.getCompoundWeightsAt(i, compounds)));
        }
        
        //Profile<DNASequence, NucleotideCompound> aligner = profile.getProfile();
        
    }
    private static void printAlignment(String name, PairwiseSequenceAligner<DNASequence, NucleotideCompound> aligner) {
        System.out.println(name);
        System.out.println(aligner.getPair());
        System.out.println("Edit distance: " + aligner.getDistance());
        System.out.println("Similarity: " + aligner.getSimilarity());
    }
}
